package matrix;

public class LUDecomp extends SSLE {

	@Override
	protected void solve(boolean buildSteps) throws Exception {
		Matrix u = a.clone();
		Matrix l = Matrix.Identity(u.width(), u.type());
		if (buildSteps) {
			steps.add("LU \\;\\; decomposition");
			steps.add("A = LU");
			steps.add(l.toLatex() + u.toLatex());
		}
		// lu decomp
		for (int i = 0; i < u.height() - 1; i++) {
			for (int j = i + 1; j < u.height(); j++)
				l.set(i, j, u.solveRows(i, j, i));
			if (buildSteps)
				steps.add(l.toLatex() + u.toLatex());
		}
		if (buildSteps)
			steps.add(l.toLatex() + u.toLatex() + "X=" + b.toLatex());
		// l(ux) = lz = b
		if (buildSteps) {
			steps.add("Backward \\;\\; substitution (1)");
			steps.add("LZ = B");
			steps.add(l.toLatex() + "Z=" + b.toLatex());
		}
		Matrix z = solveLower(l, b, buildSteps);
		if (buildSteps)
			steps.add("Z=" + z.toLatex());
		// ux = z
		if (buildSteps) {
			steps.add("Backward \\;\\; substitution (2)");
			steps.add("UX = Z");
			steps.add(u.toLatex() + "X=" + z.toLatex());
		}
		x = solveUpper(u, z, buildSteps);
		if (buildSteps) {
			steps.add("Final \\;\\; solution");
			steps.add("X=" + x.toLatex());
		}
	}

	private Matrix solveLower(Matrix l, Matrix b, boolean buildSteps) {
		Matrix z = new Matrix(1, l.height(), l.type());
		Matrix g = l.clone().augment(b);
		for (int i = 0; i < g.height(); i++) {
			z.set(0, i, g.get(g.width() - 1, i));
			if (buildSteps)
				steps.add("z_" + (i + 1) + "=" + g.get(g.width() - 1, i));
			for (int j = 0; j < i; j++) {
				z.set(0, i, z.get(0, i).doubleValue() - z.get(0, j).doubleValue() * g.get(j, i).doubleValue());
				if (buildSteps)
					steps.set(steps.size() - 1, steps.get(steps.size() - 1) + "-" + g.get(j, i).doubleValue() + "*z_" + (j + 1));
			}
			if (buildSteps)
				steps.set(steps.size() - 1, steps.get(steps.size() - 1) + "=" + z.get(0, i));
		}
		return z;
	}

	private Matrix solveUpper(Matrix u, Matrix z, boolean buildSteps) {
		Matrix t = new Matrix(1, u.height(), u.type());
		Matrix g = u.clone().augment(z);
		for (int i = g.height() - 1; i >= 0; i--) {
			if (buildSteps)
				steps.add("x_" + (i + 1) + "={" + g.get(g.width() - 1, i));
			t.set(0, i, g.get(g.width() - 1, i));
			for (int j = i + 1; j < g.height(); j++) {
				t.set(0, i, t.get(0, i).doubleValue() - t.get(0, j).doubleValue() * g.get(j, i).doubleValue());
				if (buildSteps)
					steps.set(steps.size() - 1, steps.get(steps.size() - 1) + "-" + g.get(j, i).doubleValue() + "*x_" + (j + 1));
			}
			t.set(0, i, t.get(0, i).doubleValue() / g.get(i, i).doubleValue());
			if (buildSteps)
				steps.set(steps.size() - 1, steps.get(steps.size() - 1) + "\\over" + g.get(i, i).doubleValue() + "}=" + t.get(0, i));
		}
		return t;
	}

	@Override
	public String nSteps() {
		int n = a.height();
		return (n * n) + " mul & div, " + (n * n - n) + " sub";
	}

	@Override
	public String toString() {
		return "LU Decomposition";
	}

}
